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Abstract. We discuss a recent bosonization method developed to study clean Fermi 
gases with repulsion in any dimensions. The method enables one to consider both 
density and spin excitations. It is demonstrated that due to a non-abelian structure of 
the effective theory, the spin excitations interact with each other, which leads to new 
logarithmic in temperature corrections to physical quantities. Using a renormalization 
group scheme constructed for the effective low energy field theory these logarithms are 
summed up in all orders. Temperature dependent corrections to the specific heat and 
spin susceptibility are obtained for all dimensions d = 1,2,3. 
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1. Introduction 

Landau theory of the Fermi liquid (FL) il] is by now one of the most established theories 
describing a system of interacting fermions. The main statement of this theory is that 
the low energy behavior of the system is similar to that for the ideal Fermi gas. This 
gives a possibility to discuss experimental systems omitting the interactions and using 
phenomenological constants for the effective mass of the particles, density of states at 
the Fermi surface and other physical quantities. 

A phenomenological description of the FL developed in the first works pQ was 
supported by a diagrammatic analysis |2j, which was a very good confirmation. However, 
the microscopic Landau theory of Ref. j2] is based on a very strong assumption that 
one can single out a singular particle-hole channel and sum proper ladder diagrams. 
Irreducible vertices entering the ladder diagrams should remain finite and be analytic 
in the limit of small momenta and frequencies. It is generally believed that for a Fermi 
gas with repulsion this assumption is correct and such a system should behave as the 
Fermi liquid in dimensions d > 1. 

Of course, the similarity between the FL and ideal Fermi gas cannot be exact and 
there should be corrections at finite temperatures, frequencies or momenta. Study of 
unconventional metals like high temperature superconductors, heavy fermion materials, 
etc., have revealed considerable deviations of their properties from those predicted by 
the FL. As a result, quite a few theoretical works have appeared recently where the 
Landau FL theory was discussed [31 El El El Ej ■ 

The corrections to physical quantities become especially interesting when they 
are non-analytic functions of temperature, frequency or momentum. This means 
that physical quantities like C (T) /T, where C (T) is the specific heat and T is the 
temperature, x(T), where x (T) is the spin susceptibility, etc., cannot be represented 
at low temperatures as a series in T 2 , which contrasts the ideal Fermi gas. Such non- 
analytical corrections were studied in a number of publications using diagrammatic 
expansions in the electron-electron interaction [HI El CHI [HI [121 Hal E31 CHI CEl CZl EHJ - As 
a result of this investigation, it is well known by now that in d = 3 the next-to-leading 
term in C (T) jT is proportional to T 2 lnT [HI El HOI II lj . In 2d, the non-analytical 
corrections to C (T) /T and x (T, Q) scale as T and max {T, vqQ}, respectively (Q is the 
wave vector of the external magnetic field) [131 EH 113 CHI EH EE EH] • 

The existence of the non-analytic corrections to the physical quantities is not 
accidental. In fact, all the singular corrections to the thermodynamic quantities may be 
understood in terms of contributions of low lying collective excitations, see e.g., Refs. 
[TTH IT7j . At the same time, explicit calculations with the conventional diagrammatic 
technique are not simple already in the lowest orders of the perturbation theory. 
Therefore, selecting diagrams in order to group them into the collective modes is a 
rather difficult task. 

In this paper we present a new method of calculations for a clean Fermi gas 
with a repulsion that enables us to "integrate out" electron degrees of freedom in 
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the beginning of all calculations and reduce the initial fermionic model to a model 
of low lying excitations. The method is based on using equations for quasiclassical 
Green functions and includes both the density and spin excitations. We call loosely 
our method bosonization but it differs from earlier higher dimensional bosonization 
schemes (2H1 ED EH 1221 ESI EH ESI ESI E| based on the assumption a long range electron- 
electron interaction. As a result, only density excitations were considered, while the spin 
excitations were not included in those schemes (to be more precise, the spin excitations 
are not affected by the long range interaction and therefore they were not included to 
the scheme). 

The density excitations are described by a scalar function and the effective 
interaction V c (k,u>) between them vanishes in the limit u, k — > 0. This means that the 
interaction may lead only to a renormalization of parameters characterizing physical 
quantities (unless V c (k, uj) is long ranged) but not to new effects. In other words, the 
bosonization ofj201EIliniE21ESlEllESlESlEl enables one to reduce a system with a 
long range interaction to a model of free bosons, which is very similar to writing the 
Tomonaga-Luttinger model for one dimensional electron systems (see e.g. (2jtJ)- This 
method is often referred to as bosonization. 

What we want to present now is a scheme that enables one to consider for arbitrary 
interactions both the density and spin excitations on equal footing. The relevant variable 
describing the spin excitations is a 2 x 2 matrix and our method resembles to some 
extent a non-abelian bosonization. The effective interaction between the spin excitations 
does not vanish in the limit u, k — > and leads to new logarithmic contributions to 
scattering amplitudes coming from low energies of the order of T. As a result, physical 
quantities depend on the logarithmically renormalized amplitudes, which changes their 
temperature dependence. It is important to emphasize that the "infrared" logarithmic 
divergencies that we find in the limit T —>■ exist in any dimensions (including d — 1) 
and have nothing in common with the "ultraviolet" logarithmic divergency (divergency 
originating from short distances) discussed for 2d systems long ago [2E|- The latter 
does not lead to any additional dependence on temperature and can be absorbed into 
parameters characterizing FL. 

We display in the subsequent Sections the main idea of our method and new 
results that have been obtained recently. In SecEl we make a Hubbard-Stratonovich 
transformation and derive quasiclassical equations for the density and spin excitations 
representing their solution in terms of an integral over supervectors. The interaction 
between the spin excitations gives logarithms that are summed in SecEl using a 
new renormalization group (RG) scheme. We calculate the specific heat and spin 
susceptibility in SecEJand discuss the results in Sec|SJ 

The method and the calculation of the specific heat has been presented for the 
first time in Ref. j2H]- Using this method the calculation of the spin susceptibility was 
carried out later in Ref. [SO]- Our presentation here is based on these publications but 
we concentrate rather on explaining the main steps of the derivation than on explicit 
calculations. All necessary details can be found in Refs. [23E0I- 
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2. Spin and density excitations and their contribution to thermodynamics. 

2.1. Singling out slow pairs and Hubbard- Stratonovich transformation. 

In this Section we show how one can reduce calculation of the partition function of 
the interacting fermions to calculation of the partition function of the density and spin 
excitations. It will be demonstrated that the density excitations are described by a 
model of free bosons, whereas spin excitations interact with each other. 

We start the discussion writing the original partition function Z in terms of a 
functional integral over classical anticommuting variables Xa ( x ) ( x — { r ) r } an d o~ labels 
the spin) 



Z = J exp(-S)D X D X * (1) 
The action S entering Eq. has the form 

S = So + Si n t, (2) 
where the term So, 

S = J2 J & ( x ) (—jff ~^ + ep ) x ° ^ dx (3) 

stands for the action of free fermions (e.p is the Fermi energy, m is the mass and p is 
the momentum operator). Eqs. (0121) are written in the Matsubara representation with 
the imaginary time r, such that the field variables x ( r ) r ) are antiperiodic in r 

X(r,r) = - X (r,r+1/T) (4) 

The term Si nt is Eq. (J2J) describes the fermion-fermion interaction, 

S int = \z~l J Xl( x ) Xl> 0') v(x- x) Xa> (a/) Xa ix) dxdx', (5) 

a, a 1 

where v (x — x') — V (r — r') 5 (r — r') and V (r — r') is the potential of the interaction. 

Presence of a magnetic field b acting on spin can be accounted for by adding an 
addition term Sb to the action 

S b = dx xl{x)ba aa/ Xa'(x), (6) 



Inclusion of this term is necessary for calculation of the spin susceptibility x {T)- 

The functional integral over Xa ( x ) i n Eq. (JTJ) is too complicated to be calculated 

exactly and we restrict ourselves with the case of a weak interaction. A stronger 

interaction may renormalize the coupling constants but does not seem to change the 

temperature behavior. 

In order to reduce the fermionic model to the model of the low lying excitations we 

single out in the interaction term Si nt , Eq. (jSJ), pairs x*X slowly varying in space and 

write the effective interaction Si nt as 
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Sint^l Yl J dPidP 2 dK{V 2X : (Pi) Xa (Pi + K) xl> (P 2 ) X.' (P 2 - K) 

- Vi (Pi - p 2 ) Xl (Pi) X.' (Pi + K) xl> (P 2 ) X. (P2 - K)} (7) 

In Eq. (|7j), Pi = (pi,e ni ) , where p« is the momentum and e ni = ttT (2n^ + 1) are 
Matsubara fermionic frequencies (i = 1, 2). As concerns K, it has the form K = (k,u; n ), 
where u n = 2irTn are Matsubara bosonic frequencies. 

The symbols of the integration J dPi and J dK in Eq. (J7J) read as follows 

d d p , , C .... . _ v-^ . ,. . d d k 



«(-)=TE^,(-).y^(-)=TE/«^ 



where 



/ (k) = /o (Ar ) (9) 

and k — |k|. The function / (t) has the following asymptotics: /q (t) = 1 at £ = and 
/ (t) -»• at t -»• 00. 

The function / (k) in Eq. ® is written in order to cut off large momenta k. The 
parameter r is the minimal length in the theory and we assume that r > Pp 1 . So, the 
momenta k are cut by the maximal momentum k c = r^ 1 and we avoid double counting 
when calculating the partition function Z. 

Introducing the cutoff r means that the pairs written in Eq. (J7J) vary slowly in 
space. Accordingly, we neglect the dependence of V\ and V% on the momentum k in Eq. 
(J2J). Although being smaller than the Fermi momentum pp and the Fermi energy ep, 
the cutoff k c is larger than all other momenta in the model. 

Additional decoupling in the Cooper channel is not included, since this would 
amount to overcounting of relevant scattering processes (2H]. To be short, the most 
important is parallel or antiparallel motion of the particles and one needs to consider 
only forward and backward scattering. In this limit, adding the Cooper channel would 
mean double counting. This is not so for disordered systems where all scattering angles 
are important and where one should take into account the Cooper channel |31j . 

For a short range potential we can further simplify our considerations be setting 
V2 — V(\q\ Pf). Since important momenta are close to the Fermi surface we can 
write Vi(#i2) = V(p± — p 2 ) = V(2p sin(^y-)), where #i 2 is the angle between momenta 
pi and p 2 , 6>i2 = pTpl- 

For the further development of the theory it will be crucial to separate explicitly 
interactions in the triplet and singlet channel. Making the notations 

v s (9 12 ) = v 2 - ^(0i 2 ), v t (6 l2 ) = l -v x (e l2 ) (10) 

one can represent the interaction term in a form of a sum of charge and a spin parts, 

Sint — <Sint,s + <Sint,t j (H) 

S in t,s = \ J d PidP2dq p(pi, -q)V s (9 12 )p(p 2 ,q) , (12) 
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2 . 

where the charge p(p, q) and spin densities S(p, q) are 

p{p, q) = X ] (p - |) X (p + |) , S(p, g) = x 1 " (p - |) (p + I) , (14) 

and we turned to a spinor notation x = X{)- 

In order to simplify the presentation we do not write for a while the function / 
assuming that the variables p and S are not equal to zero for small q only, which 
corresponds to a slow variation of these variables in space. 

Having written the interaction term Si n t in the form of Eq. (jllj) we next decouple 
it using a Hubbard-Stratonovich transformation with a field 

<f> a (x) = i(p a {x) + o-h n (x) (15) 

Here ip n (x) and h n (x) are real bosonic fields, so that n (r, r) = n (r, t + (3) and n is the 
direction of momentum p on the Fermi surface, n = p/|p|. The result is the following 
representation of the partition function (we omit for a while the external field b) 

Z = Jv<P W s [<p]W t [h]Z[h, cp] . (16) 

The weight functions W s , Wt are shown below in Eqs. (|2()| l2*T|) . The partition 
function Z\h, ip] describes the fermion motion for a fixed configuration of fields h, (p 

Z[<j>] = J D{ X *,x) exp(-5 e// M). (17) 
where the effective action S e ff has the form 

SefM = S + J dpdndrtxHn, r)K ( E 4^) x(r2 ' r ) eip(ri " 2) • 

Now we can write down a representation of the partition function in the presence of the 
magnetic field as a weighted integral over field configurations 



Z = J T>(p W s [(p]W t [h)Z[(j)) , (19) 
where the weights W s [<p] and W t [h] are 

W s {<p}=exp - ^ J dn 1 dh 2 dqdT(p* ni (q,T)V f r 1 (9 12 )<p m (q,T) (20) 

W t [h] = exp [ - X - J dn 1 dn 2 dqdrht li (q ) ^Vf 1 (9 12 )K 2 (q, r)] (21) 

Eqs. (|16H21|) is the final result of this subsection. We see that the original problem of 
the electron with interaction is replaced by a problem of electron motion in a potential 
and a magnetic field slowly varying in space. The condition of the slow variations 
follows from our separation into slow pairs. As we will see, at low temperature and 
weak interactions the slow variations of the fields h and tp give the main contribution 
into the physical quantities. 



Collective mode description 



7 



2.2. Quasiclassical equations 

What we should do is to calculate quantities for any h and tp and integrate these 
quantities over these fields. First, we introduce the Green functions G^, (x, x') 
corresponding to the action S e ff [0] , Eq. (fT%|). as follows 

Gt y (x, x') = Z- 1 [0] J Xa (x) xl> (x) exp (S eff [0]) D X D X * (22) 

As the fluctuating field 0, Eq. ([To]) , varies slowly in space one can derive 
quasiclassical equations for the function G^ aa , (x, x'). The method of the derivation 
is well known [32J. One should perform the Fourier transform with respect to the 
difference r — r' and assume that the Green function slowly depends on R = (r + r') /2. 
Introducing the quasiclassical Green function gt (R,t, t') in the standard way (321 ESI 

r°° d£ d 2 

gt (R,r, T J )=iJ Gl (R,r, r') £ = ^ - e F (23) 

we come to the following equation for this function 
(jt + £-iv F nV) gi(R,r,r') 

+ g* (R;r, r') n (R/) - n (R,r) ^ (R;r, r') = (24) 

where n 2 = 1, such that p F n is a vector on the Fermi surface. 

In principle, one could derive Eq. more accurately, which would produce 

additional terms containing space derivatives of the functions n (R,r) and g£ (R;r, r') 
in the second line. However, the additional derivatives would compensate infrared 
singularities we are interested in. This is the reason why we neglect them. At the 
same time, no higher derivatives arise in the first line in Eq. ()24j) and this term is exact. 

The function g& (R,r, r') must obey the antiperiodicity conditions 

gt (R,r, r') = -gt (R,r + 1/T, r') = -gt (R,r, r' + 1/T) (25) 

that follow from Eq. (@J). 

Eq. (|2*4^) is linear and therefore is not sufficient to find g^ (R,t, t') unambiguously. 
However, the same equation as Eq. (|2*4*j) can be written for the function g 2 

r l/T 

g 2 (R,r, /) = / gt (R,r, r") gt (R/, r') dr" (26) 

An obvious solution for g 2 can be written as 

g 2 (R,r, r') = cS (r - r') (27) 

where c is an arbitrary constant. It can be fixed assuming that the fermion-fermion 
interaction is present only in a finite, although macroscopic part of the space. Then, 
outside this space we come to the Green function of a free fermion gas satisfying Eq. 
(l2*7j) with c = 1. So, we come to the equation 
r i/T 

/ gt (R,r, t") gt (R,r", r') dr" = 6 (r - r') (28) 



Collective mode description 



8 



Eq. (|28|) complements Eq. (|24|) and these equations are sufficient to find the function 
(R,t, r') for any function p (R, r). After that, in order to calculate physical 
quantities, one should perform a proper averaging over p (R, r) with the weights 
W Stt [0], Eqs. (PUI l2*Tj) . In the next subsection we will show how to express the partition 
function Z[(fi], Eq. (fTTj) . in terms of the solution of these equations but now let us 
reduce Eqs. (1241 l2*Hj) to a more simple form. 

Now we come to the main point of the method proposed here. We notice that the 
quasiclassical Green function of the free fermion gas is singular and can be written as 

9on (r - r') = -iTRe [sin" 1 ttT (r - r' - id)] (29) 

where 6 — > +0. Of course, the Green function g 0n (r — r') satisfies Eqs. (J24"l 125} [2f 

As concerns arbitrary <f) n (R, r) , we look for the general solution of Eqs. 
in the following form 

g* (R,r, t') = T n (R, r) g (r - r') T~ 1 (R,r') (30) 

where T n (r,r) is a spin matrix satisfying the condition 

T n (R,r) = T n (R,r + l/T) (31) 

The representation of the Green function in the form of Eq. (j3*Uj) is a generalization 
of the Schwinger Ansatz jHl]. The form given by Eq. (|30jl . is consistent with Eq. (}28|) 
and what remains to be done is to find the proper matrix T n (R, r), such that Eq. (|2~4"j) 
is satisfied. 

A straightforward manipulation (SHj enables one to reduce Eqs. ()24| l2*%|) to the 
following form 

~ + ^ F nV R ^ M n (x) + [0 n (x) , M n (x)} = (32) 
where 

M n (x) = ^^T- 1 (x) (33) 

and the symbol [, ] stands for the commutator. 

Using the representation, Eq. (fl3j) . for the matrix n (x) and writing the matrix 
M n (x) as 

M n (x) = p n (x) + S n (x) a, (34) 

where p n (x) is a scalar function and S n (x) is a three dimensional vector, we reduce Eq. 
(JH21) to two independent equations for p n (x) and S n (x) 

~ + ^nV R ) p n (x) = -i 8 -^^- (35) 

~ + i^nVaj S n (x) + 2z [h n (x) xS n (x)] = ( 36 ) 

Eqs. dnSJ EH)) are the final quasiclassical equations that will be used for further 
calculations. We emphasize that Eqs. (JHol are obtained from Eqs. ([2H l2*H|) without 
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making any approximation. The variable p n (x) describes collective density excitations, 
whereas the variable S n (x) stands for spin ones. 

Eqs. (jHSllHEI) describing these excitations are remarkably different from each other. 
Eq. (|35|) for the density is rather simple. This is what one obtains using the high 
dimensional bosonization of Refs. [201 EU El 1221 EBl EH ESI EEl El from an eikonal 
equation. Of course, we could take into account gradients of the field (p n (x) and this 
would lead to additional terms in the L.H.S. of Eq. (|33j) . However, this does not lead 
to new physical effects. 

In contrast, Eq. (|3*Bj) is highly non-trivial due to the presence of h n (x) in the L.H.S 
of the equation. Actually, the homogeneous part of Eq. is just the equation of 

motion of a classical magnetic moment in the external magnetic field. We will see that 
the form of Eq. (}36|) will result in very non-trivial effects that will be considered later. 

The presence of the second term in the L.H.S. of Eq. (jHEJ) is a consequence of a 
non-abelian character of the variables describing the spin excitations. In this respect our 
method resembles the non-abelian bosonization well known for one-dimensional systems 

E3- 

It is important to emphasize that the variables p n (x) and S n (x) depend not 
only on the time and coordinate but also on the position of the vector n, that 
determines the position on the Fermi surface. This dependence is usual for kinetic 
equations. In the previous attempts to construct a higher dimensional bosonization 
EniiniEll221Ell23l23l2nil3 the corresponding variable arose from a "patching" of 
the Fermi surface. 



2. 3. Partition function. 



What remains to be done is to express the partition function Z [<&], Eq. (fTTT HHjl . in 
terms of the variables p n (x) and h n (x). Integrating over x, \* i n Eq. (fTTjl and using 
Eqs. (fTSl EJ) for S e ff [0] we write Z [0] in the form 

d 



Z 



cxp 



Tr I In 



P 2 



dr 2m 



+ £F + 0n (x) dx 



Eq. (}3T|) can further be represented as 



Z 



Z exp 



Tr 




(x, x\ u) n (x) dxdu 



where the Green function (r, r;r, r; u) is the solution of the equation 



d_ 



— — h ep + ucj) n (x) ) G^ (x, x'\ u) = 5 (x — x') 
2m I 



(37) 



(38) 



(39) 



We see that the calculation of the partition function reduces to the calculation of the 
Green function G^ (x, x'; u) that differs from the Green functions calculated in the 
preceding subsection by the replacement <fi (x) — * u<f)(x), where u is a parameter in 
the interval (0, 1). This means that calculating the Green function we can repeat all the 
transformations we have performed previously. 
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The Green function (r, r; r, r'\ u) at coinciding points can be written in terms 
of the quasiclassical Green function g£ as 



-iixv J gfl (r;r, t';u) dn 



(40) 



where v is the density of states on the Fermi surface (without taking into account the 
double degeneracy due to spin). 

Next, we use the representation, Eq. (|30jl . and expand the function T" 1 (r, r') in 
t' — r. The contribution from g 0n ( r _ t'), Eq. (|2l?j) . vanishes at r = r'. Using Eqs. (j33 
EH) we obtain finally 



where 



Z 



Z„ 



ZoZ^Zh 



exp 



exp 



-2iu J J p n (x, u) <p (x) dxdndu 



-2u 




S n (x, m) h n (x) dxdndu 



The functions p n (x, u) and S n (x, u) should be found from the equations 



Ju=0 (n) p n (X) 



dfn (x) 

dr 



L u (n) S n (x) = -u 



<9h n (x) 

dr 



In Eqs. 



the operator L u has the form 



L n , (n) 



d 

— — h ivp (nV r ) + 2iuh 
or 



(41) 

(42) 
(43) 

(44) 

(45) 



where the matrix h is 

/ 



h 



h- 



o 



hy \ 







(46) 



/ 



and h x , h y , and h z are the components of the vector h (ha = [h x a] for any vector a). 
The functions S n (x) and p n (x) are periodic in r 

S n (r,r) = S n (r,r + 1/T) (47) 

The accuracy of Eqs. ()20l l2~Tj) can be somewhat improved by making the 
substitution V s ,t — > F s j/i/, where f Sjt is the scattering amplitude for the singlet and 
triplet channel respectively. 

Thus, we have reduced the study of the system of the interacting fermions to 
investigation of a system of bosonic density and spin excitations. Therefore the word 
"bosonization" is most suitable for our approach. We see that the method should work 
in any dimension. At the same time, it is more general than the scheme of the high 
dimensional bosonization of Refs. [2011^10121123121112312313 because we can 
consider the spin excitations that are much more interesting than the density ones. The 
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presence of the non-trivial third term in Eqs. (J3bl l4*5)l is a consequence of the non- 
abelian character of the excitations. In contrast, the previous schemes worked only for 
abelian density excitations and it is simply impossible to write any interaction of p with 
an external field tp without using time or space derivatives. 

3. Effective field theory and renormalization group equations 

3.1. Infrared logarithmic divergences. 

Eqs. (j41H47|) . (j!9H21J) are sufficient for calculating low energy contributions to the 
thermodynamic quantities. Calculation for the density excitations is not difficult 
because the first equation (|4*4"|) can easily be solved explicitly. As concerns the spin 
excitations one may seek for the solution of the second equation (}4"4"|) expanding the 
operator L u , Eq. ()44j) . in h. It turns out that in the limit T — > terms of the 
perturbation theory for scattering amplitudes are logarithmically divergent in any 
dimension and one has to sum an infinite series. 

In order to see the origin of this divergency, let us consider the expression 
K(n,—n') = J L~ l (n;x, x') L^ 1 (— n';x', x) dx'. Using the Fourier transform in the 
coordinates and time we bring this expression in the limit T — > to the form 



If n is parallel to n' we can integrate separately over the parallel k\\ and perpendicular 
kj_ (with respect to the vector n) components of the vector k. In this case the integrand 
does not depend on kj_ and formally diverges at large |kj_|. However, we assumed that 
the momenta cannot be very large and therefore the maximal |kj_| are of the order of 
r^ 1 from Eq. (jUJ), which provides the convergence of the integral over k^. 

In contrast, the integral over uj and k\\ diverges at small values of these variables. 
These are infrared divergences and they lead to important contributions to the 
thermodynamic quantities. Estimating the value of K we come to the following 
expression 



Products of the function K (n, — n') arise when calculating scattering amplitudes with 
the help of the perturbation theory and we see that in the limit T — > one should 
sum an infinite series. This is not an easy task but the consideration simplifies if we 
reformulate the problem of calculation of the partition function in terms of a field theory. 
Then we will be able to use for summation of the logarithms a renormalization group 
technique. 

3.2. Low energy super symmetric field theory. 

Now, our task is to solve Eqs. (jSJ), substitute the solution into Eqs. (|4*Tl |4"3*|) and, 
using the obtained expression for Z [</>], calculate the integral over the fields (p n and h n . 






(49) 
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This procedure is somewhat similar to what one does in theory of disordered metals. A 
convenient way of calculations is to represent the solutions of Eqs. (jUj) in a form that 
would allow one to integrate over the fields <p n in the beginning of all the calculations. 
Integration over supervectors, Refs. |H3 EE] , is most convenient for this purpose and 
we follow now this method. As the contribution of the density fluctuations is simple, 
we consider from now on the spin excitations only. The contribution of the density 
excitations will be added in final results. 

Using supervectors ip and formulae for gaussian integration we represent the 
partition function Z h , Eq. (jlHj) . as follows [29J 



Z h = / exp[2uV2i / i) {X)T (X) dX] exp [-«S h [V 3 ]] Dip (50) 



where 



S h [V>] = -2iv J ip (X)H^j (X) dX, H = H - 2iuh (X) r 3 (51) 

The weight denominator in the gaussian integral, Eq. (JoT)|) . is absent because ip are 
supervectors. Although the general form of Eqs. (JBTfl IBTj) is simple (it is a gaussian 
integral), the detailed structure of the vectors and matrices is not. The supervectors ip 
have 48 components. The number of components comes from the necessity to consider 1) 
3 spin components (s-space), 2) bosonic and fermionic variables (g-space), 3) "particles" 
and "holes" (e/i-space) , 4) "left" and "right" motion (n-space), 5) one should double 
the number of the components to "hermitize" the space (if-space). The operator L u (n), 
Eq. ()45j) . is not hermitian and the "hermitization" of the space of the supervectors is 
necessary to provide convergence of the gaussian integral in Eq. (|50|). 

The supervectors ip are assumed to satisfy the bosonic periodicity conditions 

ip(r) = V>(r + 1/T) (52) 

We emphasize that Eq. (J52|) holds for both the bosonic and fermionic components of 
the supervectors ip. 

The generalized coordinate X contains the components 

X = (x, z) , z = (n,u) (53) 

where n is now a unit vector, n 2 = 1, parallel to the original vector n but having 
only positive x components n x > 0. Negative x-components are taken into account by 
doubling the number the components of the supervector ip (n-space). 
The operator 7i in Eq. (jHTj) can be written as 

Ho = -iv F (nV) r 3 S 3 - A^ (54) 

where Aj is the first Pauli matrix in the hermitized space (it is unity in all the other 
spaces) and the matrices r 3 and E 3 are the third Pauli matrices in the eh and n spaces 
respectively 

^0-0. ^0-°' (55) 
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The matrix h (X) has the spin structure of Eq. (|46|) (the replacement n — > nE 3 
should be made) and is unity in all other spaces. The action S^ [if)] is supersymmetric 
and is invariant under homogeneous rotations in the superspace. 

The supervector if) in Eq. (jSUJ) is conjugated with respect to if) (see for the definition 
Ref. PS])- The vector JF has the form 

F h (X) = Ti (X) + Tl (X) (56) 

where 




and 7 = x,y,z stands for the spin indices. 

We see (J5U1 l5T | [5T) ]) that the Hubbard- Stratonovich field h (X) enters both TC and T 
linearly. We write the contribution of the spin excitations Z t into the partition function 

as 

Z t = J Z h W t [h]Vh, (57) 

(c.f. Eq. (|19j)). which enables us to integrate immediately over h. 

The integration over h is gaussian and can easily be performed. As a result, we 
obtain an effective action S containing non only the free quadratic part but also cubic 
and quartic interactions 



Z t = IViP exp (-S) , S = S + S - 

J a=l,2,3 



(5f 



The interaction-independent part So equals 

So = -2iv J dX~ a (X) H ip a (X) (59) 
The three different interaction terms present in the theory can be written as 

4 r 

S 2 = -w^ X ij / dXdX 1 (60) 

i,j=l J 

x (^(x)T 3 n,<9 x jo) p^n^pri)) 

4 

S 3 = - 2\f^2iv ^ e S p~/ (61) 



x 



dXdX, (if; s (X) ursUjiPp (X)) r XjXl (Fodx^ (X x )) 
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71 



(62) 



rfXrfXx (X) nr 3 n^ 7 (X)) (X) u^H^ (X 1 )) 

Summation over spin indices a, (3, 5, 7 is implied and we use the absolutely antisymmetric 
tensor e a ^ with £123 = 1. The action S, Eqs. (I58ll62j) . is sufficient for calculation of the 
thermodynamic potential in the absence of a magnetic field. The operator T l xx , has 
the form 



rW, = 7i(nn') f{v-v')5{r-r') 



and 



7i(nni, 



72 { nrix j 




7/ 



7b 



(63) 



(64) 



(65) 



where / (r) is the cutoff function introduced in Eqs. (JH1 EJ). The operator dx in Eqs. 
has the form 



d x (a) 



1 

ud-r 



(66) 



As we will see, the most interesting contribution comes from n and ri! nearly parallel 
to each other. This justifies the notations 7/ and 7^ standing for the bare forward and 
backward scattering. 

The matrices IT equal 

n x = i, n 2 = s 3 , n 3 = Air 3 , n 4 = A 1 r 3 s 3 (67) 

and 

/ 1 1 -1 -1 \ 
1111 



A; 



1-11-1 
1-1-1 1 



(68) 



/ 



The vector Tq in Eqs. 



has the form 



V2 \ 1 



/ 



V 



-1 



(69) 



/ 



eh / H 

The vector JF projects on the bosonic sector and its presence violates the supersymmetry 
of the terms «S , a = 2, 3. 

The action S, (J59H62|) has the cubic and quartic interaction terms and looks 
quite complicated. Nevertheless, explicit calculations are not very difficult because 
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perturbation theory in the interaction is logarithmic (see Eqs. (|48( l4T?j) i n a H dimensions 
d = 1,2, 3. In order to sum up the logarithms we use a renormalization group scheme. 

At the end of this subsection we write additional terms in the action arising due 
to an external magnetic field b acting on spins. Writing these terms is necessary for 
calculation of the spin susceptibility. These terms can be brought to the form (HHI 

<Sbo = ~ v V J dx b 2 (x) (70) 
S bl = -2v4^T% r) J dX b s (x) @ Xt5 T 3 d x F ) (71) 

S b2 = 4:V £ &f3l 7] J dX bs(x) (lp X ,f3 UT 3?Px,-y) (72) 

In these expressions 

' = _itW (73) 

where the bar in V t means averaging over the full solid angle. 
3.3. Renormalization group equations and their solutions. 

We use a standard momentum shell renormalization group scheme. Separating fast and 
slow fields in the action we integrate over the fast fields and determine in this way the 
flow of coupling constants as a function of a running cutoff. In our case this amounts 
to a re-summation of the perturbation theory in the leading logarithmic approximation. 
We assume during the renormalization that the coupling constants 7 are small, 7 <C 1. 

In our case it convenient to define fast $ and slow \P fields with respect to the 
frequency only. The reason is the anisotropy in momentum. As one can see, relevant 
momenta p\\ are of the order of uj/vp, while momenta do not contribute to the 
logarithm and enter as parameters. Thus, we write 

i/>(X) = *(X) + H>(X), (74) 

where the fast fields $ have the frequencies uj in the interval, 

KLO c < \aj\ < UJ C 

while the slow ones ^ carry frequencies 

< KLU C , 

where u c is the running cut-off and k < 1. Fast modes are integrated over in the 
Gaussian approximation using averages of the form 



(...)„ = J *&(...)exp (-$,[«&]). (75) 

This results in a change in S 

5S[m\ = - In (exp (-5[¥ + $])) - S[% (76) 
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that will now be determined explicitly. 

The main object of the RG calculations that we start now is the Green function 
Go (k, n,u) corresponding to the bare action So, Eq. (j5§j) . Calculating a gaussian integral 
we obtain easily the form of the Green function in the Fourier representation 

Go (k, n,uj) = -Mv <^) = , 1 (77) 

where (...) stands for averaging with So, Eq. ()59|). 

The bare Green function, Eq. (|77jl has a non-trivial matrix structure. Therefore 
renormalized vertices are also complicated and the form of the quadratic S2 and cubic 
terms S3 in the action, Eqs. (|6UH6ip . is not most general. At the same time, the 
supersymmetric terms So and 1S4 do not change their form under the renormalization. 
Actually, So does not change in the first order and therefore we write the corresponding 
equations first for the effective vertices 7$(£) , i = 1,2,3,4, of the quartic interaction, 
where £ is the running logarithmic variable. 

A detailed calculation has been performed in Ref. R turns out that only 

the vertices 71 and 73 are non-trivially renormalized, whereas 72 and 74 remain equal 
to their bare values 7° = 7/ and 7^ = 76- At the same time, only 73 (£) enters the 
thermodynamic potential and we write here the RG equation for 73 (£) only 

= - M?)] 2 (78) 

The renormalization of the S2 and S3 terms is more complicated. In order to obtain 
a closed system of RG equations one should write a more general form of the quadratic 
and cubic interactions. This is because the renormalization of different elements in 
the particle-hole space runs in a different way. As a result there are 2 different 
cubic vertices /3 4 and (3^ and 3 quadratic vertices A 4+ , A l and A 4 " - = A^ + for all 
% = 1,2, 3, 4. The bare values of all these vertices can be found from Eqs. ()60j IHT|) and 
they are equal to 7/ for i = 1, 2 and % for i = 3,4. 

The vertices /3j and Aj do not renormalize for i = 2, 4 and remain equal to their bare 
values. Moreover, the vertices A, do not renormalize for i — 1, too. The thermodynamic 
potential is determined solely by A 3 (£) and we write down the RG equations for i = 3 
only. We remind the reader that this value of % corresponds to the backward scattering. 

The equations for the cubic vertices f3f (£) can be written as 

= -2 73 (0 /5 3 + (0 ; = -73 (0 AT (0 ; (79) 

The equation for A3" 4 " (£) takes the form 

= -2A 4 + (o 73 (0 - 2 [Pi (0] 2 ; (so) 

whereas the equation for A3 (£) is 

^K).^K)._ 2ft -« )ft+(f) . (81 ) 
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Equation for A3 ~ (£) can be written from the condition of the absence of ultraviolet 
divergencies, which follows from the invariance of the system under spin rotations. It 
takes the form 

As" (073(0= [&-(0] 2 (82) 

Eqs. (I78H82|) can easily be solved and their solutions satisfying the boundary 
conditions at £ = (when the vertices are equal to their bare values) take the form 

73 (0 = /V (0 = A3 " (0 = 77^7; (83) 

# (0 = a+- (0 = a 3 - + (0 = — fr^ ( 84 ) 

(C + 

A++ (0 = Z *> 3 5b , (85) 

where we introduced the notation 

e b = - > 0. 

7fo 

with the backscattering amplitude 7^ defined in Eq. (|B3)l. In Eqs. (|83H85|) one should 
write the final value of the £ at which the renormalization stops. Its value can be written 

as 



f (0; it, Mi; rj.) = u x u[i d f (r±) In 



1 V f 

mm I -. 



86) 



8r T 

where 9 is the angle between the vectors n and ri! and the function / is the Fourier 
transform of the function / {k\\ = 0,kj_), Eq. Q, with respect to kj_. The parameter 
/j, d equals: ^ = 2; /i 2 = 4 (p^o) _1 ; A*3 = 4?r (p F ^o) _2 • 

One can see from Eqs. (j83H85|) that all the relevant vertices decay as £ grows. This 
is usually referred to as the "zero charge" situation. If we start from small values of the 
vertices (weak coupling) they become even smaller in the process of the renormalization 
and the one loop approximation used here is sufficient for writing the final results. 



4. Thermodynamic quantities. 

4-1- Specific heat. 

Calculation of thermodynamic quantities can be performed starting with a standard 
relation for the thermodynamic potential Q (T) 

tt(T) = -TlnZ (87) 

where Z is the partition function. 

We have performed the renormalization group calculations for the case when 
the vectors n and n' of two spin excitations were close to each other (parallel or 
antiparallel motion). Only in this limit one obtains large logarithms that determine 
the renormalization of the vertices. A crucial question is whether or not this narrow 
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region of the parameters can bring an important contribution to thermodynamic or 
other physical quantities. This is not a trivial question because the system was not 
assumed to be one- or quasi-one-dimensional and one could imagine that all the effect 
of the singularities in the vertices would be washed out after the summation over the 
whole phase space. 

In fact, this almost parallel motion of the spin excitations does not contribute much 
into the thermodynamic potential fl (T) itself. Fortunately, this is not a very interesting 
quantity and what one would like to know are derivatives of the thermodynamic 
potential with respect to temperature and other parameters. One of the most important 
thermodynamic quantities is the specific heat C that can be expressed through the 
thermodynamic potential fl as 

d 2 fl 

What we need is to calculate not the thermodynamic potential fl (T) itself but the 
difference Sfl (T) , 

Sfl (T) = fl (T) — fl (0) (89) 

Using the diagrammatic method of the calculations we should be able to express the 
thermodynamic potential fl (T) in terms of sums of products of the Green functions 
over bosonic Matsubara frequencies u n 

fi(T) = ^i?K) (90) 

where R (u n ) is a function of the frequency. 

The sums of the type of Eq. ()90|) are very often divergent at high frequencies. This 
problem can be avoided calculating the quantity Sfl (T) , Eq. ([89*]) . 

Using the Poisson formula we represent Sfl (T) in the form 

mm = 2E/£*M «i>(~) < 91 > 

which improves drastically the convergence. The essential frequencies in Eq. ("9"T]) are of 
the order T and are smaller then those frequencies that form logarithms in the vertices. 

We calculate the partition function Z t using Eq. ("""*"))■ If we kept in the action 
S [ip] the super symmetric part S [ip] + 5 4 [ip], Eqs. (T59l """""]) . only, we would get zero for 
the thermodynamic potential fl. The terms <S>2 and S3 [ip] , Eqs. (JB01IHIJ), violate the 
supersymmetry and, as a result, one obtains finite contributions to fl. 

The thermodynamic potential fl (T) can be expanded in terms of the renormalized 
action 5 2 W\ and the lowest orders of the expansion take the form 

fl (T) = fl, (T) + fl 2 (T) , (92) 

where 

fl 1 (T)=T(S 2 m , r2 2 (T) = -|<( l S 2 [^) 2 > (93) 
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and (...) means averaging over ip with the Lagrangian S [ip] , Eq. (JoTJj) . The quantities 
Sfli t 2 (T) are obtained for Q 1>2 (T) by subtracting fl li2 (0) . 

It turns out that the terms 5Qi (T) and 5Q 2 (T) lead to qualitatively different 
types of contributions. In the first order in the interaction, only the forward scattering 
contributes. As we have found in the previous subsection, the part of S 2 corresponding 
to the forward scattering is not renormalized and one may use just Eq. (joTIj) for it. 
Then, we come to the expression 

where Ao differs from r , Eq. ©, by a numerical coefficient. Actually, Ao should be 
of the order of the wave length if one wants to make estimates for the initial model 
of the interacting fermions. The factor 3 in the numerator is due to the fact that the 
spin-1 excitations have 3 projections. A similar contribution comes also from the density 
excitations but, of course, with the factor 1. Eq. describes a contribution of the 
interaction to the coefficient in front of the linear dependence of the specific heat on 
temperature. 

The part 5Q 2 consists of the part further renormalizing the coefficient C/T of the 
linear dependence on T and a part giving corrections non- analytic in T 2 . A general 
expression for 5Q 2 (T) can be brought to the form 

5n 2 (T) = -6 mn £; I ^ exp (-i^ -r,\^ 

x jdn x dn 2 J -^ 3 Y(9;k ± ,k ll )P d (u;Xn l ,n 2 ) (95) 

where 9 is the angle between the vectors ri! and n 2 . The main interesting contribution 
will come from small 8, which justifies the decomposition of the momentum k into 
perpendicular k^ and parallel kn with respect to 111,2 components. 
The function Y kj_, fcii) defined as 



u\u 2 du\du 2 (9 



r(0;ki,A;||) = f [ 
Jo Jo 

x I [A^ _ (0; ui, u 2 \ k ± , fcy)] 2 + (0; ux, u 2 \ k ± , k\\) A^ - (9; u x , u 2 \ k ± , k\\] 

is the most important entry in the final expression for the specific heat. The vertices 
A 3 should be taken from Eqs. (j84H86j) . 
The formfactor 

^ / \ jiuj + v F n 2 ) (iu-vpnx) 

V d a;, ni, n 2 = — — ■ (97) 

[luj — Vpn- 2 ) [iuj + Vpn-i) 

depends on the dimensionality of the system and it describes basically the free 
propagation of the two spin excitations in almost opposite directions. The non-analytic 
contributions originate from the small region of the phase space |ni — n 2 | <C 1. 

The result of the calculation depends on the dimensionality of the system but the 
non-trivial corrections exist in both d = 2 and d = 3. The details of the calculations 
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can be found in Ref. j2H] but the final results for the non-analytic corrections 5C to the 
specific heat look for d = 2 and d = 3 as 

^o s { w|!in ¥ + ^r$ iLM - 2), i ( "» 

where A" (T) = fi d 1b In (s /T) (e = vp/r ~ e F and Li 2 (x) = J2T=i xk /k 2 is the 
polylogarithm function. The first terms in the circular brackets in Eqs. (J951 l9"9"j) describe 
the density excitations and 7^ is the coupling constant for the backward scattering of 
these excitations. The second terms originate from the spin excitations. The final result, 
Eqs. (jnEEU), was written for a special choice of the cutoff function f (r±), Eq. (fHBJ) 

/(r ± ) = ^ 1 exp(-r ± /r ) (100) 

where Qd is the d — 1 dimensional solid angle (Q\ = 2, f2 2 = 27r). 

In the limit of not very low temperatures when X (T) <C 1, Eqs. (|HH| EHJ) take a 
simpler form 



3vr 4 / T \ 3 



[7, P ] 2 + 3 7b 2 ) (101) 



^ = -^(^m(^)( [7fl V«) (102) 

Eqs. (I10H I102|) agree with results obtained by a direct diagrammatic expansions (see 
for a recent discussion Refs. [T7J[TH]. At the same time, the general equations (J9"2~l l9*3~j) 
have been obtained for the first time in Ref. [21] using the bosonization scheme. We 
see that in the limit T — > the contribution of the spin excitations vanishes and only 
the density excitations contribute to the specific heat. 

The bosonization method presented here gives a possibility to make calculations 
also for one-dimensional systems. Surprisingly, the calculation of the specific heat Id 
is more complicated than for d = 2, 3. This is because the function Vd (u, n%, n 2 ), Eq. 
(|97|) . is exactly equal to unity and the contribution of Eq. to the specific heat 

vanishes. In order to obtain a non-vanishing contribution one should consider terms of 
higher orders in the effective vertices. As a result, one comes to the following expression 
for the correction to the specific heat in d = 1 

rji -1 

5C d =i = — -3 (103) 

vf (l + 2 7b ln^) 3 

The correction to the specific heat for Id can be extracted from the exact solution for 
spin chains of Ref. jSZj- This correction agrees with our result, Eq. (|1U3|) . and we 
see that our supersymmetric low energy theory reproduces all the previously known 
physical effects despite the fact that the intermediate degrees of freedom differ from the 
conventional bosonization. 
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4-2. Spin susceptibility 

Calculation of the spin susceptibility has been performed in a recent work jHOj- The 
bosonization scheme of Ref. [2H] was used and it was shown that the logarithmic 
contributions arise as well. An external magnetic field b leads to the additional terms 
Sbo,Sbi and Sb2, Eqs. fJ70H72j) . in the effective action. 

We proceed as before integrating over the fast variables $, Eq. (J71|) and thus 
deriving renormalization group equations. In principle, one should renormalize not 
only the terms S 2 ,S 3 ,S4, Eqs. (j60H62|) . as it was done previously, but also the terms 
Sbo,Sbi and Sb2, Eqs. (I70H72|) and S , Eq. (JoTlj) . Again, the calculations should be done 
separately for d = 1 and d = 2,3. 

It turns out that in the higher dimensions d = 2, 3 the terms Sbo, Sbi and Sb2, and 
Sq do not change under the renormalization and what remains to be done is to express 
the susceptibility x m terms of the renormalized vertices 7(£), (3 (£) and A(£). The 
computation is somewhat more cumbersome than the one for the specific heat and there 
are many contributions that have to be combined together. As for the specific heat, only 
the vertices with % = 3 contribute and one can express the susceptibility through 73 (£), 
(3 3 (0 and A 3 (0- 

Using the solutions of the RG equations, Eqs. (|83H85|) . one can bring the expression 
for the susceptibility to a form containing the function Y [9; k±, k\\) , Eq. (|9fij) . Choosing 
again the cutoff function / in the form of Eq. (jlOU|) we bring the temperature dependent 
correction b~\ (T) to the spin susceptibility to the form [HO] 

fe(^8«l&!^™ (104) 

WT) = Q^. (ff& ™-£r>™' (105) 

where rj is given by Eq. (|73*|l . The susceptibility xf ^ s ^ ne Pauli susceptibility in d 
dimensions. 

In the limit X (T) 1 one can neglect the logarithmic contributions and obtain 
the following form of the corrections 5xd 

5x^(T) = 2r?>fi—j& (106) 

Ef 

^3(T) = y«(0 X i°2 3 (107) 

The linear dependence on T of the correction x (T) , Eq. ()106|) , agrees with those 
obtained by conventional methods ^Sl El HS1 QS|- At the same time, the proportionality 
of Sx (T) to T 2 in d = 3 is analytical in T 2 and Eq. ()107|) describes a renormalization 
of a coefficient in front of the T 2 . As the T 2 -term is present already in the temperature 
dependence of the susceptibility of the ideal Fermi gas, the correction is not very 
interesting. The first non- analytical term T 2 ln(e^/T) in the temperature dependence 
in d = 3 is proportional to 7^. 
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Calculation of the susceptibility for one dimensional systems is somewhat more 
involved because on should consider corrections to the terms Sbo,Sbi and Sb2, Eqs. (I7U1 - 
l?2*|) and So, Eq. (|59j) . Nevertheless, one can proceed in a rather straightforward way 
and, as a results, the following dependence of the susceptibility is obtained [3U| 



This result agrees the with the one obtained long ago j3H] by a completely different RG 
method developed for the electron problem in Id. So, Eq. (jl08|) serves as one more 
check of the bosonization method reviewed here. 

5. Discussion. 

In the previous sections we reviewed the new method of bosonization for a clean Fermi 
gas in any dimensions suggested in Ref. [2H| and further used in Ref. [SHI- In contrast 
to previous attempts [23 I2D 13 121 123 121 123 123 13 we do not restrict ourselves with 
the case of a long range electron-electron interaction and include into the scheme spin 
degrees of freedom. This enables us to consider not only density excitations but also 
the spin ones. In contrast to the density excitations that can be described by a scalar, 
the spin excitations are described by 3-component vectors and one may speak about a 
new version of the non-abelian bosonization. 

The non-abelian character of the effective theory leads to the non-trivial interaction 
between the spin modes. Making the perturbation theory in this interaction new 
logarithmic contributions (diverging in the limit T — > 0) to vertices were discovered 
and summed up using a RG scheme. As a result of this consideration we have 
found temperature dependent corrections SC (T) /T to the specific heat and to the 
susceptibility 5x(T) in all dimensions, Eqs. (|§H1 E3 11031 11041 I105J) . The results for 
d = 1 agree with those obtained earlier by completely different methods. In higher 
dimensions, d = 2, 3 the lowest order of the expansion of our formulae in the coupling 
constant agree with known results obtained previously using conventional diagrammatic 
techniques, see e.g. [El H3 H3 EJII3 EH) • All these agreements serve as a good check 
of our approach. 

The new contributions to the specific heat and susceptibility originate from an 
almost parallel motion of the spin excitations. Although both the forward and 
backward scattering amplitudes are renormalized, only the backward scattering enters 
the thermodynamic quantities. 

In the language of conventional diagrams the logarithmic contributions to the 
thermodynamic quantities come from both Cooper and particle-hole loops because they 
originate from quasi-one-dimensional processes. The fact that the forward scattering 
amplitude drops out from the final results corresponds to what happens in Id. 

Although we have agreement in the limiting cases with almost all the results we 
could compare with, there is a disagreement with an old work j^H] , where an instability 
of the Fermi liquid against superconducting pairing with a high angular momentum was 



S X d=i (T) 



(108) 
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found. This is especially strange because the main contribution to the formation of the 
superconductivity comes from almost parallel electron motion such that only forward 
and backward scattering enter the superconducting critical temperature. This is just the 
region that we considered. As we do not see any such effect, this can mean that either 1) 
the accuracy of our method is not sufficient (we summed the terms like (7 In (ef/T)) 12 , 
whereas (j 2 In (ep/T)) n were summed in Ref. jSH]) or 2) summing the Cooper ladder 
as was done in Ref. [SHj is a bad approximation for the quasi-one-dimensional process. 
We believe that 2) is the explanation of the contradiction but a more careful analysis 
should be performed. 

We have developed and applied the bosonization scheme for the simplest model of 
the interacting Fermi gas. One can add other terms to the Hamiltonian of this model 
to take into account different features relevant for experimental systems. We hope that 
our scheme may be useful in this study. 
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